library(rstan)

cat("\nCCES 2008 PID\n\n")

# Set to the local directory of the replication data
setwd("~/Downloads/Replication_Data/")
n_cores <- 8
options(mc.cores = n_cores)

L2008D <- readRDS("Data/CCES-L2008D.rds")
L2008R <- readRDS("Data/CCES-L2008R.rds")

##### Stack data to feed to Stan

model_data_2008D <- list(n_items = length(unique(L2008D$j)),
                        n_respondents = length(unique(L2008D$k)),
                        n_obs = nrow(L2008D),
                        k = L2008D$k, j = L2008D$j, y = L2008D$y)

model_data_2008R <- list(n_items = length(unique(L2008R$j)),
                        n_respondents = length(unique(L2008R$k)),
                        n_obs = nrow(L2008R),
                        k = L2008R$k, j = L2008R$j, y = L2008R$y)

##### Compile to C

model_am <- stan_model(file = "AM_Sigma.stan")

##### Sample from posterior

exclude_pars <- c("beta_raw", "alpha_raw", "zeta_raw")

# Democrats
n_items <- length(unique(L2008D$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2008D <- sampling(model_am, data = model_data_2008D,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2008D, "Fitted_Models/posterior_cces_2008D.rds", compress = "xz")


# Republicans
n_items <- length(unique(L2008R$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2008R <- sampling(model_am, data = model_data_2008R,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2008R, "Fitted_Models/posterior_cces_2008R.rds", compress = "xz")

